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ABSTRACT 

Context. Recent detections of disks around young high-mass stars support the idea of massive star formation through accretion rather 
than coalescence, but the detailed kinematics in the equatorial region of the disk candidates is not well known, which limits our 
understanding of the accretion process. 

Aims. This paper explores the kinematics of the gas around a young massive star with millimeter-wave interferometry to improve our 
understanding of the formation of massive stars though accretion. 

Methods. We use Plateau de Bure interferometric images to probe the environment of the nearby (~1 kpc) and luminous (~20000 
Lq) high-mass (10-16 Mq) young star AFGL 2591-VLA3 in continuum and in lines of HDO, Hi^^O and SO2 in the 115 and 230 
GHz bands. Radiative transfer calculations are employed to investigate the kinematics of the source. 

Results. At -0.5'' (500 AU) resolution, the line images clearly resolve the velocity field of the central compact source (diameter 
of ~ 800 AU) and show linear velocity gradients in the northeast-southwest direction. Judging from the disk-outflow geometry, the 
observed velocity gradient results from rotation and radial expansion in the equatorial region of VLA3. Radiative transfer calculations 
suggest that the velocity field is consistent with sub-Keplerian rotation plus Hubble-law like expansion. The line profiles of the 
observed molecules suggest a layered structure, with HDO emission arising from the disk mid-plane, H2^^0 from the warm mid- 
layer, and SO2 from the upper disk. 

Conclusions. We propose AFGL 2591-VLA3 as a new massive disk candidate, with peculiar kinematics. The rotation of this disk 
is sub-Keplerian, probably due to magnetic braking, while the stellar wind may be responsible for the expansion of the disk. The 
expansion motion may also be an indirect evidence of disk accretion in the very inner region because of the conservation of angular 
momentum. The sub-Keplerian rotation discovered in our work suggests that AFGL 2591-VLA3 may be a special case linking 
transition of velocity field of massive disks from pure Keplerian rotation to solid-body rotation though definitely more new detections 
of circumstellar disks around high-mass YSOs are required to examine this hypothesis. Our results support the idea that early B-type 
stars could be formed with a circumstellar disk from the point of view of the disk-outflow geometry, though the accretion processes 
in the disk need to be further investigated. 

Key words. Stars: massive - Stars: formation - ISM: individual objects: AFGL 2591 - Accretion, accretion disks - ISM: kinematics 
and dynamics 



1. Introduction 

Massive stars play an important role in cycling and mixing ma- 
terial between stars and the interstellar medium in galaxies. 
However, our understanding of how high-mass stars (M > 8 
Mq) form is still far from complete (Zinnecker & Yorke 2007). 
Observationally, it is challenging to gain information from high- 
mass young stellar objects (YSOs) because high-mass stars usu- 
ally form in groups at large distances (typically a few kpc), re- 
quiring high angular resolution to resolve individual sources. 
High-mass YSOs are embedded in regions with high extinctions 
(~ 10-100 Ay) for ~ 15% of their lifetime, where complex pro- 
cesses, such as gravitational interaction, powerful outflows and 
stellar winds, high accretion rates and ionizing radiation fields 
make the interpretation of observations difl&cult. Moreover, high- 
mass YSOs evolve very fast (~ 10^ yrs) and reach the main se- 
quence while still in the embedded phase. Structures such as cir- 
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cumstellar disks have an even shorter lifetime and are harder to 
detect. 

Accretion via a circumstellar disk is suggested to be the most 
eflfective way to overcome radiative pressure and build up a high- 
mass star (e.g., Krumholz et al. 2009). The detection of massive 
bipolar outflows with high mass loss rates of ~ 10""^ Mq yr"^ 
around young high-mass stars implies the presence of circum- 
stellar disks with comparable mass accretion rates (e.g., Henning 
et al. 2000; Beuther et al. 2002; Zhang et al. 2005). To date, 
only a dozen or so B-type young stars have been proposed to 
have circumstellar disks (e.g., Cesaroni et al. 1997, 1999, 2006; 
Shepherd & Kurtz 1999; Beuther et al. 2005), with disk masses 
less than or comparable to the masses of the host stars. In con- 
trast, around young 0-type stars, only large rotating toroids are 
found with masses greater than stellar masses (Beltran et al. 
2005;Furuya et al. 2008). 

The presence of circumstellar disks around low-mass young 
stars is firmly established, and our current observational and the- 
oretical understanding is reviewed by Williams & Cieza (2011) 
and Armitage (2011), respectively. In contrast, the applicability 
of a scaled-up version of low-mass star formation via disk accre- 
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Table 1. Log of observations 
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03 Feb 2006^ 


19 Feb 2007 


13 Mar 2007 


15 Mar 2007 


Wavelength 


3 mm + 1.3 mm 


3 mm 


1.3 mm 


1.3 mm 


On- source time (hours) 


5.3 


4.5 


5.0 


4.0 


Configuration 


6Aq 


6Aq 


6Bq 


6Bq 


Bandpass caHbrator 


3C454.3 


3C273 


3C345 


3C273 


Gain calibrator 


2005+403, 2013+370 


2005+403, 2013+370 


2005+403, 2013+370 


2005+403, 2013+370 


Flux calibrator 


MWC349 


MWC349 


MWC349 


MWC349 



Notes. 

Antenna 4 missing due to tuning problems. 



tion for massive stars has not been firmly established (Zinnecker 
& Yorke 2007), although promising cases have been reported 
(e.g., Keto & Zhang 2010). Clearly, more observational work is 
needed to establish the kinematics of material associated with 
high-mass YSOs. 

In this paper, we present observations of massive disk can- 
didate, AFGL 2591-VLA3, a young early B-type star in the 
Cygnus X region. The distance to this object is uncertain, with 
reported values ranging from 0.5 to 2 kpc (0.5-2 kpc as dis- 
cussed in van der Tak et al. 1999) and recent measurements of 
water masers even suggesting 3.3 kpc (Rygl et al. 2011). In this 
work, we assume a distance of 1 kpc in order to aid compari- 
son with previous publications and discuss how our conclusions 
change for a distance of 0.5 or 3 kpc. At 1 kpc, the total lu- 
minosity of AFGL 2591 is ~ 2 x 10"^ Lq (Lada et al. 1984). 
Compact radio emission is associated with the source (named as 
VLAl, VLA2 and VLA3; Campbell 1984; Trinidad et al. 2003). 
Modeling of its free-free emission implies the stellar mass is 
about 16 M© (van der Tak & Menten 2005). An envelope mass 
of only 40 M© within 30000 AU around the source is derived 
through comprehensive radiative transfer modeling of molecular 
lines by van der Tak et al. (1999). 

The geometry of the system can be derived from bipolar 
^^CO outflows at arcnunute scale, which extend in east (red)- 
west (blue) direction (position angle (P.A.) -270° -325°) ob- 
served at resolutions of 14'' to 21'' by Mitchell et al. (1992) 
and Hasegawa & Mitchell (1995). Observations of ^^CO at res- 
olution of 14" suggests that the P.A. of the outflow is about 
280° (van der Tak et al. 1999), confirmed by i^-band bispec- 
trum speckle imaging of AFGL 2591 at a resolution of ~ 0.2" 
showing multiple loops on the western side of the source with 
P.A. ~ 260° (Preibisch et al. 2003). Precession of the outflow 
axis is also implied. Dense gas in the outflow walls traced by 
CS and HCN shows a consistent morphology (Bruderer et al. 
2009). Based on this orientation of the outflow, a north- south 
(P.A. ~ 0°) orientation of the equatorial plane of AFGL 2591- 
VLA3 is inferred. 

Observations of HDO li,o-li,i, H2^^0 3i,3-22,o and SO2 
12o,i2-lli,ii with the Plateau de Bure interferometer (PdBI) of 
the Institut de Radio Astronomic Millimetrique (IRAM)^ re- 
veal a velocity gradient in the north(east)-south(west) direction 
across a barely resolved source with a diameter of ~ 800 AU, 
suggestive of a circumstellar disk (van der Tak et al. 2006). 
However, the kinematics could not be studied due to the limited 
angular resolution (~ 1"). Here we present new PdBI observa- 
tions toward AFGL 2591 at a higher resolution of ~ 0.5" (~ 500 
AU at 1 kpc), which resolve the kinematics. This paper is orga- 
nized as follows. Section 2 summarizes the PdBI observations 
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and data reduction. Section 3 presents the results. In Sect. 4, we 
describe a simple model of the kinematics. Discussions and con- 
clusions are given in Sect. 5 and Sect. 6, respectively. 



2. Observations 

The source AFGL 2591 was observed with PdBI in 2006-2007 
in 4 tracks. The 3 mm receiver was tuned to cover HDO li^o- 
li,i at 80578.3 MHz (hereafter HDO) and the 1.3 mm receiver 
was tuned to cover H2^^0 3i,3-22,o at 203407.5 MHz (hereafter 
H2^^0) and SO2 12o,i2-l li,ii at 203391.6 MHz (hereafter SO2). 
Table 1 presents a log of the observations and Table 2 summa- 
rizes the observed lines. The combined projected baselines range 
from 26 k.1 (-10") to 206 k.1 in length at 3 mm and from 48 k.1 
(~5") to 525 k/1 in length at 1.3 mm. The phase tracking center 
was chosen at Qf(2000) = 20^29"^24^.87, ^(2000) = 40°11'19'.'50. 
The nominal LSR velocity adopted for observations was -5.5 
km s"^ The quasars 2005-h403 and 2013-h370, which are about 
5° away from AFGL 2591, were frequently observed for phase 
and gain calibrations in all four observations. 3C273, 3C345 and 
3C454.3 were used as the bandpass calibrators, and MWC349 
with a total flux density of -1.3 Jy at 203.4 GHz during the 
time of the observations, as flux calibrator. The flux error is es- 
timated to be within 20%. The spectral resolutions are 39.0625 
kHz per channel or 0.145 km s~^ per channel at 3 mm band, and 
78.125 kHz per channel or 0.115 km s~^ per channel at 1.3 mm 
band. Data reduction was conducted at the IRAM headquarters 
in Grenoble, using the GILDAS^ software. 



Table 2. Observed molecular lines 









E ^ 


Beam Size^ 


Molecule 


Transition 


(MHz) 


(K) 


rx-°) 


HDO 


1 1,0-1 1,1 


80578.3 


47 


1.23 x0.72; 26 


H2^^0 


3 1,3-22,0 


203407.5 


204 


0.51 X 0.33; 22 


SO2 


12o,i2-lli,ii 


203391.6 


70 


0.51 X 0.33; 22 



Notes. 

^""^ Rest frequency taken from JPL molecular spectroscopy 
(http://spec.jpLnasa.gov/). 

Energy of upper level. 

Synthesized beam using natural weighting. 
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80.6 GHz 



203.4 GHz 

"^n ^ r- 

VLA 3-NE 




RA (J2000) 

Fig. 1. Natural- weighted continuum images in the region of 
AFGL 2591. (Left) 80.6-GHz continuum. Contours are at 3, 5, 
10, 20, 30... cr with 1-cr of 0.15 mJy beam"^ Sources within the 
field following Campbell (1984) and Trinidad et al. (2003) are 
marked on the map. (Right) 203.4-GHz continuum. Contours are 
at 3, 5, 10, 30, 50... a with l-a of 0.4 mJy beam'^. 
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Fig. 2. Visibility analysis of the continuum at 203.4 GHz. (a) and 
(c) show the observed visibilities (filled circles) and model vis- 
ibilities (grey area). The observed data are vector- averaged for 
the purpose of better visualization. Only the averaged structure 
in circular symmetry can be observed in the plot. Full model vis- 
ibilities are plotted in order to show the non-circular symmetric 
structure of the source, (b) and (d) represent a graphical view 
of the visibility fits of VLA3 shown in (a) and (c), respectively. 
The red plus sign is the peak position of the emission. Blue el- 
lipse and filled circle stand for the Gaussian component and the 
point source. 



3. Results 

3.1. Continuum emission 

The natural-weighted continuum images of AFGL 2591-VLA3 
at 80.6 GHz and 203.4 GHz are presented in Fig. 1. The angular 
resolutions are 1724 x a.'72, RA. 28° at 80.6 GHz and a.'51 x 
a'33, RA. 18° at 203.4 GHz. With an arcsecond beam, VLA3 
is unresolved at 80.6 GHz. In contrast, at 203.4 GHz VLA3 is 



well resolved with a subarcsecond beam and shows a diamond 
shaped morphology in the north-south direction. In addition to 
VLA3, VLAl and VLA2 are detected within the field of view 
at 80.6 GHz. At 203.4 GHz, VLAl is likely resolved out due to 
missing short spacings. A previously unknown source (VLA3- 
NE) located at -2.5'' north-eastern to VLA3 is detected at 5-cr 
at 203.4 GHz. 



3.1 .1 . Measurements of positions, sizes and total flux 
densities 

The positions, sizes and total flux densities of the sources in the 
region of AFGL 2591 are derived in the visibility domain and 
summarized in Table 3. Judging from the images in Fig. 1, all 
sources are fit with point sources, except VLAl at 80.6 GHz 
and VLA3 at 203.4 GHz, where Gaussians are used. The vector- 
averaged visibility plot in Fig. 2 shows a compact component 
on baselines longer than 200 kA for VLA3 at 203.4 GHz. We 
therefore fit VLA3 with a 20 mJy point source plus a 65 mJy, 
1.0'' X 0.6" (~ 800 AU diameter at 1 kpc) Gaussian source. 
Including the flux densities of VLA2 and VLA3-NE, the total 
flux density within the field of view at 203.4 GHz is about 93 
mJy which is consistent with the visibility measurement (Fig. 
2). 

Compared to van der Tak et al. (2006), our 80.6-GHz obser- 
vations measure 100%, 12% and 45% of the total flux densities 
of VLAl, VLA2 and VLA3, respectively, while the new 203.4- 
GHz observations detect 0%, 7% and 44% of the total flux den- 
sities of VLAl, VLA2 and VLA3, respectively. Considering the 
difl'erent array configurations (limited by the antenna shadowing 
eff'ect, the shortest baselines of the observations made by van der 
Tak et al. (2006) are 8 k/1 at 3 mm and 20 k.1 at 1.3 mm) and the 
absolute uncertainty of the flux density, the new data are consis- 
tent with optically thin free-free emission, optically thick dust 
or free-free emission and optically thin dust emission for VLAl, 
VLA2 and VLA3, respectively, as discussed by Trinidad et al. 
(2003) and van der Tak et al. (1999, 2006). For the newly de- 
tected source VLA3-NE, more sensitive high-resolution obser- 
vations at difl'erent frequencies are needed to derive its nature. 

3.1 .2. Estimation of densities and masses 

For the "dust" source VLA3, we derive the beam averaged H2 
column density from the measured peak intensity of 30 mJy per 
beam at 203.4 GHz. The H2 column density at the emission peak 
can be estimated as 



A^(H2) = 



2mii^},KyBy(T^) ' 



(1) 



where 5'y is the peak flux density, a is the gas-to-dust ratio 
(100), Qb is the beam solid angle, mu is the mass of atomic hy- 
drogen, Ky is the dust opacity per unit mass and By(T^) is the 
Planck function at dust temperature T^. We apply the value of 
/Cy(1.3mm) 1.0 cm^ g"^ as suggested by Ossenkopf & Henning 
(1994) for gas densities of 10^ - 10^ cm"^ and coagulated dust 
particles with thin ice mantles. The beam averaged H2 column 
densities for from 100 to 300 K (see Sect. 3.2.4 for the valid- 
ity of this assumption) are 0.5 - 1.7 x 10^^ cm"^. The gas mass is 
estimated from the integrated continuum flux density of 85 mJy 
via 
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gas 



Fycfa 

KyByiT^) ' 



(2) 
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Table 3. Characteristics of Continuum Emission 





Aa^ 


A^« 


n b n b 
'^major i^minor 


P.A.^ 




F^(vdT06)^ 


Source 


n 


n 


n n 


n 


(mJy) 


(mJy) 


80.6 GHz 


VLAl 


-3.67 


-4.49 


2.28 2.13 


-18 + 7 


62 + 1 


61 + 1 


VLA2 


-3.64 


1.31 






1.1 + 0.1 


9± 1 


VLA3 


0.11 


-0.04 






7.2 + 0.1 


16 + 2 


203.4 GHz 


VLAl 












61 + 2 


VLA2 


-2.75 


0.94 






4.7 + 0.4 


65 + 4 


VLA3^ 


0.03 


0.00 


0.96 0.60 


6 + 2 


65 + 2 


194+1 




0.16 


-0.05 






20+1 




VLA3-NE 


1.80 


1.82 






4.7 + 0.4 





Notes. Positions, source sizes and total flux densities are derived in visibility domain. 

Position offset with respect to q'(2000) = 20^29^^24^.87 and ^(2000) = 40°iri9'.'50. 
ib) pwHM source size and position angle. 

Total flux density. The absolute flux uncertainty is ~ 20%. 
^"^^ Total flux density measured by van der Tak et al. (2006) in the visibility domain. 

Model with a Gaussian plus a point source. See Fig. 2 and text for more discussions. 



Table 4. H2 column density, mass and number density of VLA3 



Dust Temperature 





100 K 


200 K 


300 K 


Note 


A^(H2) / W cmr' 


1.7 


0.8 


0.5 




Mg^J Mq 


0.34 


0.16 


0.11 


d - \ kpc 


n(U2) 1 10^ cm-3 


22.3 


10.8 


7.2 


if size = 800 AU 


Mgas/ Mq 


3.03 


1.48 


0.98 


d - 3 kpc 


w(H2) / 10^ cm-3 


7.5 


3.6 


2.4 


if size = 2400 AU 



where Fy is the total flux density of dust emission and d is the 
distance to the source. At a distance of 1 kpc and temperatures 
of 100-300 K, the mass is 0.1-0.3 M©, significantly less than 
the stellar mass of VLA3 (;^10-16 M©, Lada et al. 1984; van 
der Tak & Menten 2005). If the distance is 3 kpc, the gas mass 
is modified to 1-3 M©, still less than the scaled stellar mass 
of > 20 M©. Regardless of the distance, the inferred mass of 
the surrounding gas around VLA3 is much less than the stellar 
mass. With the further assumption of a spherical source of 0.8'' 
in diameter, as measured from the visibilities of VLA3 at 203.4 
GHz, H2 number densities of 7 - 22 x lO'^ cm"^ are derived for 
a distance of 1 kpc, and of 2 - 7 X 10^ cm~^ for a distance of 
3 kpc. These densities are lower limits if the source structure is 
flattened. The inferred gas number density is consistent with the 
envelope density profile within radii of 400 AU (10^~^ cm"^) as 
derived by van der Tak et al. (1999). We sunmiarize these results 
in Table 4. 



3.2. Line emission 

Compact emission of HDO, H2^^0 and SO2 is found toward 
VLA3 (Fig. 3(a)-(c)). With the improved angular resolutions of 
our new observations (natural weighting; 3mm: r.'23x0'.'72, P. A. 
26°; 1.3mm: a'51 x a.'33, P.A. 22°), all three molecular images 
are resolved and show a nearly round shape, unlike the diamond 
shape with a centrally peaked feature seen in the 203.4-GHz con- 
tinuum. At Icr noise level, no significant line emission is found 
toward other sources in AFGL 2591. 



Table 5. Source positions and size parameters of VLA3 derived 
from molecular species. 





Aof^ 


A(5^ 


n b 
C'major 


n b 
Cminor 


PA.^ 


Vgrad P.A.^ 


Molecule 


n 


n 


n 


n 


n 


n 


HDO 


0.13 


-0.08 


1.03 


0.86 


48 + 15 


54 + 12 


H2^^0^ 


0.07 


0.05 


0.84 


0.39 


15 + 7 


40 + 7 




0.14 


0.06 










S02^ 


0.09 


-0.08 


0.98 


0.92 


-2 + 5 


33 + 9 




0.13 


0.12 











Notes. Positions and source sizes are derived in visibility domain. 

Position off^set with respect to a(2000) = 20^29"^24^.87 and ^(2000) = 
40nri9'.'50. 

(b) pwHM source size and position angle. 

Position angle of the velocity gradient measured from blue- and red- 
shifted components. See also Fig. 3. 

Model with a Gaussian plus a point source. 



3.2.1 . Measurements of positions and sizes 

We derive the positions and the sizes of the line emission in 
the visibility domain (Fig. 4). In the vector- averaged visibility 
plots, the amplitude of the visibility decreases up to 200 k/l and 
is more or less flat onward. This feature is similar to what we ob- 
served in the continuum visibility at 203.4 GHz, which implies a 
point source surrounded by an extended source. As a result, we 
adopted an elliptical Gaussian plus a point source in the visibility 
fit of H2^^0 and SO2. For HDO, we used a Gaussian only in the 
fit since our baselines at 80.6 GHz only extend to 200 k/l. The 
source sizes derived in this approach are 0.4''-1.0''. The lower 
end of the result is estimated from H2^^0, which could be due to 
the large dispersion of the visibilities or reflect a smaller source 
size traced by this higher excitation line {E^ = 204 K; 47 and 
70 K for HDO and SO2, respectively). The size measurements 
of the molecular lines are consistent with those derived from the 
continuum, indicating that the dust and the molecular gas have a 
similar spatial distribution. The locations of the emission peaks 
diflfer by much less than the synthesized beam and are negligi- 
ble. We summarize the results of the visibility analysis in Table 
5. 
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moment I moment I moment I 

-1 0.1 <Vlsr< 0.6 km s ' -5.35 <Vlsr< 0.6 km s ' -1 0.1 <Vlsr< -5.35 km s ' 




Fig. 3. (a)-(c) Zero moment maps of HDO, H2 O and SO2 around VLA3, respectively. The contours are at -3, 3, 6, 9...cr for all 
three maps, where Icr noise levels are 10, 37 and 65 mJy km s"^ for HDO, H2^^0 and SO2, respectively. The green crosses mark the 
positions of emission peaks. The green ellipses represent FWHM sizes derived from visibility fit. The green filled circles indicate 
the positions of the additional point sources in the visibility fit of H2^^0 and SO2. The natural-weighted clean beam sizes, shown in 
blue filled ellipses, of HDO and H2^^0/S02 are 1723 x a.'72, P.A. 26° and a.'51 x a.'33, P.A. 22°, respectively, (d)-(l) First moment 
maps of HDO, H2^^0 and SO2 around VLA3, respectively. For all nine maps, the color scales are set to be the same. Dashed lines 
represent the isovelocity contours with an interval of 0.25 km s"^ (m)-(o) Molecular emission integrated over the velocity ranges 
of -5.35 < Vlsr < 0.6 km s"^ (red contours) and -10.1 < Vlsr < -5.35 km s"^ (blue contours) for HDO, H2^^0 and SO2, 
respectively. Contours are plotted with 10% interval of peak value of each velocity component. The inferred position angles of the 
velocity gradients across VLA3 are plotted as green dashed lines. 



3.2.2. Velocity gradients across VLA3 

The new dataset allows us to spatially resolve the velocity distri- 
bution across VLA3. In Fig. 3 (d)-(f), velocity gradients in the 
northeast-southwest direction can be seen with the red-shifted 
part toward the southwest side and the blue-shifted part toward 
the northeast side. Apart from the diff'erence in angular resolu- 
tion, HDO and H2^^0 show a similar distribution in isovelocity 
contours. Detailed investigation of H2^^0 image reveals that the 
isovelocity contours near the systemic velocity (green part) ex- 
hibits a mirrored S shape, while SO2 shows a normal S-shaped 
distribution of isovelocity contours. The diff'erence implies that 
HD0/H2^^0 and SO2 trace diff'erent kinematics of VLA3, where 
SO2 may be dominated by the east- west outflow. First moment 
maps plotted in the velocity ranges of -5.35 < Vlsr < 0.6 km 
s-i (Fig. 3(g)-(i)) and -10.1 < Vlsr < -5.35 km s'^ (Fig. 30)- 
(1)) support our speculation that SO2 emission is dominated by 
the outflow since the extreme velocities originate near the source 
center, while HDO and H2^^0 are less afl'ected by the outflow 
and presumably dominated by the equatorial motion. We derive 
the position angles of the velocity gradients across the source 
seen in Fig. 3 (d)-(f) by integrating over the velocity ranges of 
-5.35 < Vlsr < 0.6 km s'^ and -10.1 < Vlsr < -5.35 km 
s"^ for each line and overplotting them in Fig. 3 (m)-(o), respec- 



tively. The position angles of velocity gradients are calculated 
from the emission peaks of the blue- and red-shifted components 
derived from Gaussian fit. For HDO, H2^^0 and SO2, the posi- 
tion angles are 54 ± 12°, 40 ± 7° and 33 ± 9°, respectively (Table 
5). 

The orientation of the observed velocity gradients (northeast- 
southwest) is not perpendicular to the large scale outflow (east- 
west) (Mitchell et al. 1992; Hasegawa & Mitchell 1995; van der 
Tak et al. 1999) with a difl'erence of ~ 45° in position angle. We 
rule out the possibility of contamination from the bipolar outflow 
since the red-shifted emission of the molecular emission lies to 
the southwest, which is inconsistent with the red lobe of the 
outflow in the east. High-resolution ^-band imaging (Preibisch 
et al. 2003) shows that the orientation of the large-scale outflow 
remains unchanged on small scales, ruling out contamination by 
a small-scale outflow to explain the observed velocity pattern. 
A more likely interpretation for the observed velocity pattern is 
rotation in the equatorial plane around AFGL 2591-VLA3. For 
purely azimuthal motion we would expect the velocity gradi- 
ent to be perpendicular to the bipolar outflow direction. Section 
4 further explores the eff'ect of radial motion in the equatorial 
plane that afl'ects the orientation of the velocity gradient. 

To further quantify the observed velocity gradients, we plot 
the position-velocity maps with the position angles defined in 
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Fig. 4. Observed visibilities of HDO, H2^^0 and SO2 (filled cir- 
cles) overplot with model visibilities (grey area). The observed 
data are vector-averaged corresponding to the circularly sym- 
metric structure. Full model visibilities are plotted in order to 
show the non-circular symmetric structure of the source. 



Fig. 3 (m)-(o) (Fig. 5). Clearly, in all cases a linear velocity gra- 
dient is found, which is very diff'erent from the typical Keplerian 
rotation found in the protoplanetary disks around low-mass stars 
(e.g., Simon et al. 2000). We derive (projected) velocity gradi- 
ents of 5.5 ± 0.4, 10.2 ± 1.3 and 12.4 ± 0.5 km s"^ arcsec"^ or 
km s-i per 1000 AU (at 1 kpc) for HDO, H2^^0 and SO2, re- 
spectively, suggesting that SO2 is more aff'ected by the east- west 
outflow while HDO is less afl'ected. To minimize the contami- 
nation of the outflow and highlight the motion in the equatorial 
plane of VLA3, we also derive the velocity gradients by taking 
the position- velocity cut with P. A. = 0° . A consistent strength of 
linear velocity gradient (~ 15 km s"^ arcsec"^) is derived from 
all three molecules, indicating that the efl'ect of outflow is indeed 
minimized. 



3.2.3. Molecular line profile 

Spectra of HDO, H2^^( 
sion peak (Aa = -0.15'', A^ = 0.02'') in a single synthesized 
beam (1'.'23 x a.'72, P.A. 26°) are shown in Fig. 6. A triple- 
peaked line profile is observed for HDO, while the H2^^0 and 
SO2 line profiles are smooth Gaussians with weak hints of sec- 
ondary peaks coincident with those seen in HDO. We use three 



Gaussians to decompose the HDO line profile (Table 6) and find 
that Gaussians with FWHM of ~ 1.20 ± 0.05 km s"^ at Vlsr 
of -6.89 ± 0.03, -5.35 ± 0.34 and -3.92 ± 0.26 km s'^ best 
reproduce the observed spectrum, although significant residuals 
remain near -8.0 and -3.0 km s"^ The velocity diff'erences be- 
tween the blue-shifted and systemic components (1.54 + 0.34 km 
s"^), and, the red-shifted and systemic components (1.43 ± 0.43 
km s"^) are comparable, which suggest that the red- and blue- 
shifted components are likely from the same structure undergo- 
ing systematic motion, not just a spatial coincidence of three dif- 
ferent gas clumps. To decompose the H2^^0 and SO2 line pro- 
files, we fit three Gaussians fixed at Vlsr of -6.89, -5.35 and 
-3.92 km s"^ derived from HDO line profile (Table 6). For an 
FWHM line width of 1.45 km s"^ this fit reproduces most of the 
line profile, with only significant residuals at extreme velocities 
(near -2.5 and -8.5 km s"^ for H2^^0, and near -2.0 and -9.0 
km s"^ for SO2). Based on the line profiles and the observed 
velocity gradient across the source, we suggest that the velocity 
structure of VLA3 consists of three components: two exhibiting 
motions (azimuthal and radial) in the equatorial plane and one 
static component. 



3.2.4. Molecular excitation 

We estimate the excitation conditions of VLA3 based on the 
line ratios of the chemically related molecules HDO and H2^^0. 
If we assume optically thin emission and local thermodynamic 
equilibrium (LTE), the total column density at a given excitation 
temperature can be derived via 



Mr 



2X)4Waot^^ 20 -2 
X 10 cm ^, 



(3) 



where W is the integrated intensity in Jy beam"^ km s"^ ^« and 
6t are the FWHM widths of the synthesized beam in arcsec, Qtot 
is the partition function, is the upper energy level in K, Tex 
is the excitation temperature in K, and gK are the spin and K 
degeneracies, respectively, S is the line strength, ji is the dipole 
moment in Debye, and y is the rest frequency in GHz. The source 
filling factor is assumed be unity and background emission is 
neglected. In other words, the excitation temperature can be ex- 
pressed as a (nonlinear) function of the column density ratio of 
HDO and H2^^0, if HDO and H2^^0 share the same excitation 
temperature. Once the excitation temperature is estimated, the 
column density can be derived via Equation (3). 

The column density ratio A^tot[H2^^0]/A^tot[HDO] can be ex- 
pressed as 



Mot[H2 ^'O] 



Mot[H20] 



Mot[H2 ^'O] 



A^tot[HDO] Mot[HDO] Mot[H20] 



(4) 



For the column density ratio of H2O and HDO, we adopt the 
abundance ratio [H20]/[HD0] of 1400-2000 for the inner re- 
gion of the envelope around VLA3 (assuming this ratio can 
be applied to the observed rotating structure), which is derived 
based on the "jump" model of ID radiative transfer calculations 
(see Table 11 of van der Tak et al. 2006). We assume the col- 
umn density ratio of H2^^0 and H2O to be 0.002 based on the 
[^^0]/[^^0] ratio (Wilson & Rood 1994). Therefore, the column 
density ratio A^tot[H2^^0]/Mot[HDO] is 2.8-4.0. The total col- 
umn densities of HDO and H2^^0 as a function of Tq^ are shown 
in the upper three panels of Fig. 7, while the column density ra- 
tio A^tot[H2^^0]/A^tot[HDO] as a function of T^^ can be seen in 
the bottom three panels of Fig. 7. From the column density ra- 
tio plots, a range of excitation temperatures for a given velocity 
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Fig. 5. (a)-(c) Position- Velocity maps of HDO (P.A. 54°; see Fig. 3 (m)-(o)), H2^^0 (PA. 40°) and SO2 (PA. 33°). Linear velocity 
gradients are fitted to the maps and shown as dotted lines, (d)-(f) Same as (a)-(c) except the position angle of the position- velocity 
cut for all three molecules is 0° . 
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Fig. 6. Molecular spectra taken at the mean emission peak after 
convolving all data to the same angular resolution of r.'23x(y.'72, 
P.A. 26° . Fitted Gaussians are overplotted in color. Residuals are 
shown in dots. 



component can be estimated by finding the intersections of the 
derived curves of the column density ratio from the new observa- 
tions (in purple) and the horizontal lines of the adopted column 
density ratios from Equation (4) (in green). We summarize the 
molecular excitation condition of VLA3 in Table 7. 

We estimate the excitation temperatures for the blue- shifted, 
systemic and red-shifted components to be 120-165 K, 155-240 
K and 120-165 K, respectively, or 90-290 K, 110-540 K and 
90-290 K, respectively, if the error in absolute flux calibration 



is considered. We find that the blue- and red-shifted components 
have similar excitation temperatures, implying again that both 
components are from a single source structure such as the ob- 
served equatorial gas component around VLA3. There is also an 
indication that the systemic component might be warmer. The 
overall column densities counting all three velocity components 
(without absolute error in flux density) for HDO, H2^^0 and SO2 



are 1 x 10^^, 4 x 10^^ and 1 x 10^^ cm"^, respectively. The frac- 
tional abundances are estimated to be 5 x 10"^, 2 x 10"^ and 
6 X 10"^ for HDO, H2^^0 and SO2, respectively, adopting the 
203.4-GHz continuum peak flux at the resolution of 1^23 xa.'72, 
P.A. 26° and using = 200 K in the calculation of Equation 1 . 

We note that the line ratio analysis of HDO and H2^^0 
discussed above is sensitive to the assumed abundance ra- 
tios of HDO/H2O and i^O/^^O. For example, if the assumed 
A^tot[H2^^0]/A^tot[HDO] is smaller (larger) by a factor of two, the 
corresponding temperatures would be less than 100 K (greater 
than 250 K) for all three velocity components. Alternatively, if 
we adopt a temperature of 100 K (the ice evaporation temper- 
ature), the derived column density ratio A^tot[H2^^0]/A^tot[HDO] 
would be less than 7, consistent with van der Tak et al. (2006). 
Only observations of multiple transitions of HDO and H2^^0 can 
further constrain the excitation conditions of the inner region of 
VLA3. 

Alternatively, the excitation conditions on large scales near 
VLA3 can be estimated from the HDO single-dish data pub- 
lished by van der Tak et al. (2006) with the population dia- 
gram method (Goldsmith & Langer 1999). Including the ef- 
fects of beam dilution and line opacity, the excitation temper- 
ature, total column density and source size are solved self- 
consistently. The observed data (red open circles) and best fit 
model (green crosses) are plotted in Fig. 8. The excitation tem- 
perature is estimated to be 130 ± 30 K with a column density of 
~ (1.0 ± 0.1) X 10^"^ cm"^ distributed uniformly in a source with 
angular size of 17 ± 2''. Line opacities of all the transitions used 
are estimated to be less than 0.006. The excitation temperatures 
derived from the single-dish data are consistent with the ones de- 
rived from our new interferometric data. The total HDO column 
density including all three velocity components is about 3x10^^ 
cm"^ if we observe the same structure (roughly 0.8'' - 1.0'') seen 
by the interferometer with a 17" beam which is consistent with 
the value derived from single-dish data if missing flux is con- 
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Fig. 7. H2^^0 and HDO column density and column density ratio plots. The left, middle and right columns represent the plots for 
the blue- shifted, systemic and red- shifted components, respectively. In the first row, total column densities at diff'erent excitation 
temperatures are plotted in blue (H2^^0) and red (HDO) thick curves. Their column density density ratios, shown in the second 
row, are plotted in thick purple curves. The thin curves represent the errors. In addition, the H2^^0 over HDO column density ratio, 
derived from van der Tak et al. (2006) and Wilson & Rood (1994), are plotted as green horizontal lines. The excitation temperatures 
for a given component can be read out by finding the intersection of the green lines and thick purple lines. The thick vertical dashed 
lines mark the derived temperatures for ratio 2.8 and 4.0. Here we derive the excitation temperature without considering the error of 
absolute flux density. Derived temperatures and column densities are summarized in Table 7. 

Table 6. Line parameters of HDO, H2^^0 and SO2 derived from fitting of three Gaussian components. 
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km s~0 






km s~^) 






kms-i) 


HDO 


-6.89(3) 


1.20 (5) 


0.072(3) 


-5.35(34) 


1.20(5) 


0.080(3) 


-3.92(26) 


1.20(5) 


0.094(3) 




-6.89 


1.45 


0.544(19) 


-5.35 


1.45 


0.807(19) 


-3.92 


1.45 


0.712(19) 


SO2 


-6.89 


1.45 


1.64(11) 


-5.35 


1.45 


1.74(12) 


-3.92 


1.45 


2.20(11) 



Notes. Errors in units of the last decimal are given in the parentheses. Values without parentheses are fixed in the fitting. 
Line center velocity. 
FWHM line width. 
Integrated intensity. 



sidered. As a result, both single-dish and interferometric data 
imply that HDO is present in both the outer envelope and the in- 
ner dense region toward VLA3 with a temperature greater than 
-100 K. 

4. Kinematics of the equatorial region of VLA3 

4. t. Evidence of rotation and outward radial motions 

The new (sub)arcsecond resolution PdBI observations resolve 
the velocity field in the inner region of AFGL 2591-VLA3. 
Clear linear velocity gradients in the northeast- southwest direc- 
tion are observed in HDO, H2^^0 and SO2 lines. Comparison 
of the distributions of isovelocity contours of these lines reveals 
that HDO and H2^^0 best reflect the kinematics of the equato- 
rial region of VLA3, while SO2 is more afl'ected by the east- west 
outflow. In addition, a triple-peaked line profile is observed in 
HDO. Given the bipolar outflow geometry in the east- west direc- 



tion (Mitchell et al. 1992; Hasegawa & Mitchell 1995; Preibisch 
et al. 2003), which defines a north-south equatorial plane, we 
propose that the observed lines trace rotation in the inner en- 
velope around VLA3. We suggest that radial motions alter the 
velocity field's position angle from the expected north-south di- 
rection to the observed northeast- southwest direction. Consider 
an inclined disk-like structure at P. A. 0°, with rotation and radial 
motion but without vertical motion (as expected for a model in 
vertical hydrostatic equilibrium), four possible velocity patterns 
exist (as shown schematically in Fig. 9), depending on whether 
the rotation appears clockwise or counterclockwise as seen by 
the observer, and whether the radial motion are inward or out- 
ward. Only the combination of counterclockwise rotation and 
outward radial motions reproduce the observed velocity pattern 
(c.f. Fig. 3 and 9). 
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Table 7. Excitation temperatures and column densities. 
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Notes. 

Excitation temperatures estimated from column density ratios. See Sect. 3.2.4 for details. 
Total column density derived under LTE and optically thin assumptions. 

With respect to H2. Assuming a dust temperature of 200 K, Mot[H2] of 2.35 x 10^^ cm~^ is derived from the 203.4-GHz continuum at a resolution 
of ir23 X a'72, P.A. 26°. The uncertainty in Mot[H2] is not considered. 
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Fig. 8. HDO rotation diagram plotted with the data published 
by van der Tak et al. (2006). Red open circles are the data ob- 
served with JCMT and IRAM 30m. The green crosses represent 
the best-fit model from population diagram analysis (Goldsmith 
& Langer 1999). 
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4.2. A toy model 

To explore if a kinematical model as suggested in the previous 
section can reproduce the observed velocity pattern and molecu- 
lar line spectra, we constructed a toy model and calculate the re- 
sulting emission using the axisymmetric radiative transfer code 
RATRAN (Hogerheijde & van der Tak 2000). Our model includes 
a flattened "disk-like" structure surrounded by a static spheri- 
cal envelope and a bipolar outflow cavity (Fig. 10). We adopt a 
parametrized description of the density, temperature, and veloc- 
ity field, and assume LTE and optically thin conditions for the 
line emission. Given the critical densities of the observed transi- 
tions (7 X 10^, 3 X 10^ and 5 x 10^ cm"^ for HDO, H2^^0 and 
SO2 at 200 K, respectively), and the densities derived in Sect. 
3.1.2, the excitation is likely in LTE. We scale the model inten- 
sity to match the observed line strengths, which means that we 
can derive conclusions about the velocity field but not the den- 
sity or temperature. Our only aim is to show that our simple toy 
model can reproduce the observed velocity pattern and spectral 
line shapes. We do not perform any global model optimization 
in density, temperature and velocity profiles. Instead, for fixed 
sets of density and temperature profiles, we conduct a simplified 
optimization to the parametrized velocity field. Our derived pa- 
rameters are therefore not best-fit solutions but rather indicative 
values. 







Fig. 9. Schematic view of the velocity fields combining radial 
motion and rotation. 
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Fig. 10. Gas density (left) and gas temperature (right) distribu- 
tions of the source in the radiative transfer models, including a 
disk-like structure, a spherical turbulent envelope, and an out- 
flow cavity. 



4.3. Details of the model 

The density profiles for the disk-like structure (for simplic- 
ity, we use "disk" hereafter in this section) and the envelope 
are assumed to follow power-law distributions, respectively, as 
ndisk(r,0) = An^oir/r^y sin^ (6) and ^env(^) = n^oir/r^y, 
where Hqo is the density at radius (the outer radius of the disk). 
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Fig. 11. Model moment 1 maps of the combined velocity profiles 
in azimuthal and radial directions. The model position- velocity 
maps are included in each panel (bottom-left: P. A. 0°, bottom- 
right: RA. 45°). In all cases, linear velocity gradient can be 
inferred. However, only Keplerian-like rotation (~ r~^-^) plus 
Hubble-law like radial expansion (~ r^^) can reproduce the 
observed mirrored S-shape isovelocity contours seen in H2^^0 
(c.f., Fig. 3(e)). 
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Fig. 12. Simplified parameter optimization of our kinematic 
model. The grey scales represent the similarity of model and ob- 
served HDO triple-peaked line profile (lower values for better 
match). The solid curves with values in unit of km s"^ arcsec"^ 
indicate the measured strength of linear velocity gradient at P. A. 
0° in the grid models. The open circle highlights our best match- 
ing parameter set (vrq = 10 km s"^ arcsec"\ ^ = 0.35). 



A is the disk-to-envelope density ratio at rd in the mid-plane, 6 
is the angle measured from the polar axis, / is the flattening pa- 
rameter and a is the power-law index. In the disk density profile, 
regions with densities lower than An^o are set to have zero den- 
sity, resulting in a dumpling-like, flattened structure. The overall 
density in the model is the sum of ^disk and ^env For the tem- 
perature profiles in the disk and the envelope, we assume power- 
law distributions, respectively, as rdisk(^) = BT^oiR/r^)'^ and 
Tcrwir) = TQo(r/r^)~^ (i.e., we adopt a vertically isothermal 



model for the disk) where R denotes the radial direction in cylin- 
drical coordinates, Tqq is the temperature at r^, B is the disk-to- 
envelope temperature ratio, and JS is the power-law index. The 
gas temperature for a given position is the weighted mean of the 
disk and envelope values: T = (^disk7^disk+^env7^env)/(^disk+^env). 
The distance to the source is assumed to be 1 kpc. The inner ra- 
dius of the model is set to 12.5 AU (half the size of the HII re- 
gion; van der Tak & Menten 2005) and the outer radius to 400 
AU (= Td, based on the size derived from our 203.4 GHz contin- 
uum). In addition, a bipolar outflow cone with an opening angle 
of 60° is included in the model by setting the density inside this 
cone to be zero. We describe the velocity field in the disk by 
three orthogonal components as radial (vr), vertical (v^) and az- 
imuthal (v^). The velocity field in the static envelope is set to be 
purely turbulent, inspired by the observed spectra showing emis- 
sion peaks near systemic velocity. The gas velocity for a given 
position is the weighted mean of the disk and envelope values: 
^ ~ ('^env^env + ^diskVdisk)/(^disk + ^env). The inclination of the 
source is set to be 30° (van der Tak et al. 2006) with the blue- 
shifted outflow cone pointed to the west. Model parameters are 
chosen to ensure that the LTE and optically thin assumptions are 
valid. We adopt A = 10, ^eo = 1 x 10^ cm'^ a = 1.5, / = 5, 
B = 1, Teo = 100 K, = 1, a molecular fractional abundance 
of 1 X 10"^^ and a turbulent line width of 0.5 km s"^. A de- 
tailed description of the parametrized velocity field is given in 
the next section. We scale the line intensity to the observed val- 
ues, and only analyze the velocity signatures (moment 1 maps 
and position- velocity maps) and spectral line profiles. We note 
that this scaling process makes the exact profiles of density and 
temperature irrelevant. 

4.4. Evidence of sub-Keplerian rotation and Hubble-law like 
radial expansion traced by HDO and H2^^0 

The analysis of the position- velocity maps of the observed lines 
toward VLA3 reveals that the velocity field is characterized by a 
linear velocity gradient (Fig. 5), which suggests that the equa- 
torial region is undergoing solid-body rotation. However, the 
analysis demonstrated in Fig. 9 implies that an extra outward ra- 
dial component is required to reproduce the observed northeast- 
southwest velocity pattern (Fig. 3 (d)-(f)). Therefore, a direct 
link between the observed linear velocity gradient and the in- 
ferred solid-body rotation may not be trivial. In our toy model, 
we consider two types of velocity profiles in the azimuthal and 
radial directions for the disk-like structure: one is proportional 
to r"^-^ and the other follows r^^ Although all the combinations 
of the velocity field can reproduce the triple-peaked line signa- 
ture seen in HDO and the linear velocity gradient seen in the 
position- velocity maps, only Keplerian-like rotation (v« ~ ^-0-5) 
plus Hubble-law like radial expansion (vr ~ r^^) can reproduce 
the mirrored S-shape isovelocity contours observed in H2^^0, 
which is less aff'ected by the outflow than SO2 (Fig. 11; c.f.. Fig. 
3(e)). An additional velocity component in the z direction may 
redistribute the isovelocity contours to match the observed SO2 
moment 1 map. In our toy model, the outflow component is not 
included due to the complexity of outflow morphology and ve- 
locity field. Therefore, we use HDO and H2^^0 data to constrain 
the kinematics of the equatorial region of VLA3. To summarize, 
we adopt the parametrized velocity for the disk-like structure as: 

Vr = VRor 

Va^^^GMJr (5) 
v, = 0. 
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Fig. 13. (Left) Best matching model spectrum (a), moment 1 maps (b)-(d) and position- velocity map (f) for HDO. The color scales 
and contour levels for moment 1 maps are identical to Fig. 3 (d)-(l). The dashed line in the position-velocity map represents a linear 
velocity gradient with strength of 15 km s"^ arcsec"^ (Right) Best matching model for H2^^0. 



where vrq is an analog of the Hubble constant, ^ is a constant 
between and 1, G is the gravitational constant and M^ = 16Mq 
(van der Tak & Menten 2005) is the stellar mass. To estimate 
the strength of the two free parameters vrq and ^ in our kine- 
matic model, we perform a simplified parameter optimization 
which utilizes the peak velocities of the HDO triple-peaked line 
profile (-6.7, -5.4 and -3.9 km s"^ for the blue-shifted, sys- 
temic and red- shifted components, respectively), and the mea- 
sured linear velocity gradient of ~ 15 km s"^ arcsec"^ at P. A 
0° as the constraints. We constructed a grid of models for HDO 
with vro between 4 and 20 km s"^ per 1000 AU and ^ between 
0.2 and 0.5. For each set of parameters (vrq, ^), we compute A^ = 
^(observed peak velocities - model peak velocities)^, which is 
smaller for a better match. Figure 12 shows that the models with 
the lowest A^ values span a diagonal in (vrq, ^)-space, indicat- 
ing that the model is degenerate. To break the degeneracy, we 
compute the strength of the linear velocity gradient at P.A. 0° 
in our grid model and plot it as the solid curves in Fig. 12. By 
matching the observed value, we conclude that based on our toy 
model, the azimuthal motion of the disk-like structure traced by 
HDO can be characterized by sub-Keplerian rotation (about 1/3 
the strength of Keplerian rotation with a stellar mass of 16 M©) 
and the radial motion is the accelerated expansion with the ana- 
logue Hubble constant of 10 km s"^ per 1000 AU at 1 kpc. The 
best matching model spectrum, moment 1 maps and position- 
velocity map for HDO can be found in Fig. 13 (a)-(e). 

To model the H2^^0 data, we use the best matching solution 
for HDO (vro = 10 km s"^ arcsec"\ ^ = 0.35) as the starting 
point. Certainly, the turbulent line width of 0.5 km s"^ is too nar- 
row for H2^^0 since its line profile is Gaussian-like with hints 
of secondary velocity components. We therefore adopt a turbu- 
lent line width of 1.0 km s"^ in the model. The best matching 
model for H2^^0 can be seen in Fig. 13 (f)-(j). Though not per- 
fectly matching the observations, our toy model of kinematics 
is able to reproduce both HDO and H2^^0 data qualitatively. In 
our toy model, we found that the kinematics of the proposed 
disk-like structure around VLA3 is best described by two ve- 
locity components: sub-Keplerian rotation and Hubble-law like 
expansion. We also found that the conclusion of rigid-body ro- 
tation across the source judging directly from the linear velocity 



gradient observed in position-velocity maps may not be made if 
the observed velocity gradient does not align well with the equa- 
torial plane defined by the outflow. 

5. Discussion 

5.1. An embedded disk-like structure in AFGL 2591-VLA3? 

Our (sub)arcsecond PdBI images of HDO, H2^^0 and SO2 re- 
veal a clear velocity gradient in the northeast- southwest direc- 
tion across the source AFGL 2591-VLA3 (Sect. 3.2.2). Together 
with the shape of the HDO line profile (Sect. 3.2.3), and the 
symmetric temperature distribution of the blue- and red- velocity 
components (Sect. 3.2.4), we propose that the source AFGL 
2591-VLA3 is surrounded by a disk-like structure undergoing 
rotation and expansion. The kinematics of the equatorial region 
is further characterized as sub-Keplerian rotation plus Hubble- 
law like expansion as shown by radiative transfer modeling of 
a toy model (Sect. 4). In the models for HDO, a microturbulent 
line width of ~ 0.5 km s"^ matches the triple-peaked line profile, 
which appears small for the violent environment of massive star- 
forming regions (Zinnecker & Yorke 2007). On the other hand, 
a microturbulent line width of 1.0 km s"^ is needed to match 
the almost single-peaked line profiles seen in H2^^0. We ex- 
pect that an even larger turbulent line width is required to model 
SO2, judging from the observed SO2 spectrum. This suggests 
that these molecules trace diff'erent regions within the disk-like 
structure. 

We hypothesize that the disk-like structure has a layered dis- 
tribution with HDO close to the mid-plane, H2^^0 in the warm 
mid-layer and SO2 in the upper disk layers (Fig. 14). In ad- 
dition, significant emission from all three molecules originates 
in a static turbulent envelope surrounding the disk-like struc- 
ture, especially for H2^^0 and SO2 where it dominates the line 
profiles (Table 7). In the disk-like structure, we find tempera- 
tures of about 120-165 K, well above the evaporation temper- 
ature of water ice (Fraser et al. 2001). The HDO molecules 
may be formed on grain surfaces (Tielens 1983; Parise et al. 
2005), and released into the gas phase when the disk-like struc- 
ture is heated by shocks, stellar radiation and winds. Compared 
to HDO, H2^^0 originates from a more turbulent region at larger 
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Fig. 14. A cartoon showing the idea of a layered disk-like struc- 
ture around AFGL 2591-VLA3. HDO is close to the cooler, qui- 
escent mid-plane. H2^^0 is distributed in the warm mid-layer 
and SO2 originates in the upper disk layers. 



height with increased turbulence in the disk- like structure. In this 
picture, the presence of H2^^0 but absence of HDO in the mid- 
layer of the disk-like structure could be due to increased local 
water formation because of warm gas chemistry as opposed to 
evaporation of deuterated ices formed at low temperature. SO2 
is a good shock/outflow tracer in massive star-forming regions 
(Schilke et al. 1997). Judging from the facts that the observed 
SO2 line profile has significant emission from line wings, and 
that a relatively large velocity gradient is measured in SO2 com- 
pared to HD0/H2^^0, we suggest that SO2 traces a part of the 
disk-like structure where the outflow interacts with the upper 
layers. The concentric velocity patterns seen both in the blue- 
and red-shifted part of moment 1 maps of SO2 (Fig. 3(i) and (1)) 
likely support the idea since the most extreme velocities occur in 
the region close to the star (unlike HDO and H2^^0) as expected 
for an outflow origin. 

5.2. Kinematics of the proposed disk-like structure and 
comparison to other disk sources 

Our toy model suggests that the proposed disk-like structure 
around VLA3 undergoes sub-Keplerian rotation which cannot 
be directly inferred from the observed position- velocity maps. 
We speculate that a magnetic field slows the rotation down be- 
low pure Keplerian, which may help to transport the angular 
momentum outward in VLA3. Similar to the formation of low- 
mass stars, the magnetic field has been suggested to play an im- 
portant role in the formation of high-mass stars as well (Girart 
et al. 2009; Beuther et al. 2010). In strongly magnetized clouds 
during collapse, magnetic braking is thought to regulate the for- 
mation of disks around stars (e.g., Basu & Mouschovias 1994; 
AUen et al. 2003; Meflon & Li 2008). Certainly in our case, 
the magnetic field does not fully dominate the kinematics of 
VLA3, otherwise a solid-body like velocity field would be ex- 
pected in the inner region of the disk-like structure (c.f. Galli 
et al. 2006, Fig. 2). From OH maser observations, Hutawarakorn 
& Cohen (2005) measure a maser disk with radius about 750 
AU and the magnetic field strength of a few mG toward AFGL 
2591-VLA3, which may be sufficient to slow down the rotation. 
Subarcsecond dust polarization measurements of the magnetic 



field toward VLA3 and hence deriving the mass-to-flux ratio 
would be helpful to test our hypothesis if the magnetic field in- 
deed regulates the kinematics of the equatorial regions of young 
high-mass stars. 

In addition to sub-Keplerian rotation, we found that Hubble- 
law like expansion occurs in the equatorial region of VLA3. 
The overall kinematics is similar to the SiO maser disk around 
Orion-IRc2 in which Keplerian rotation and decelerated expan- 
sion characterize the kinematics (Plambeck et al. 1990). The ex- 
panding motions in the disk could be due to the interaction be- 
tween the stellar radiation and winds with the disk-like structure, 
providing the forces to push the gas outward (e.g., Hollenbach 
et al. 1994; Yorke & Welz 1996; Richling & Yorke 1997). On 
the other hand, the presence of an expanding velocity compo- 
nent in the disk-like structure could be the result of accretion 
through the mid-plane of the very inner region so that the ex- 
pansion helps to conserve the angular momentum in the equa- 
torial region. The infrared speckle imaging of AFGL 2591 by 
Preibisch et al. (2003) suggests that the outflow or wind is vari- 
able in time thus the accretion in the disk may not be contin- 
uous but episodic. In addition, little free-free emission is de- 
tected toward the source, indicating an evolutionary stage before 
the development of an ultracompact HII region (van der Tak & 
Menten 2005). As a result, the main accretion phase might have 
ceased although models show that massive young stars can con- 
tinue to gain mass through the disk even after an ultracompact 
HII region has been developed (e.g., Keto 2007). Future high- 
resolution observations of multiple HDO lines and other dense 
gas tracers, such as N2D^ and CH3CN, would be useful to con- 
strain the physical condition and the kinematics of the disk-like 
structure around VLA3 better and understand its evolutionary 
stage of high-mass star formation. 

Although massive disks showing pure Keplerian rotation 
traced by various molecules are reported such as the cases in 
IRAS 20126-h4104 (C^'^S; Cesaroni et al. 2005), AFGL 490 
(C^^O; Schreyer et al. 2006) and NGC 7538S (H^^CN; SandeU 
et al. 2003), cases showing non-Keplerian rotation were ob- 
served as wen such as AFGL 2591 (HDO, H2^^0; this work) 
and IRAS 18089-1732 (NH3; Beuther & Walsh 2008). In ad- 
dition, massive disk candidates have a wide range of mass and 
size distributions derived via various tracers as summarized by 
Cesaroni et al. (2007). Despite the fact that the detections of 
massive disks to date are still limited, Keplerian rotation in the 
disks is typically found toward young massive stars with masses 
about 10 Mq such as the cases quoted above. For stellar masses 
about 16 Mq, massive disks start to show non-Keplerian mo- 
tion such as rigid-body rotation reported in the literature. For 
even higher stellar masses, the rotating structures (toroids) typ- 
ically show solid-body rotation and are gravitationally unstable 
(Beltran et al. 2005; Furuya et al. 2008). The sub-Keplerian ro- 
tation discovered in our work suggests that AFGL 2591-VLA3 
may be a special case linking transition of velocity field of mas- 
sive disks from pure Keplerian rotation to solid-body rotation 
though definitely more new detections of circumstellar disks 
around high-mass YSOs are required to examine this hypothesis. 

5.3. Effect of the uncertainty in distance 

The distance to AFGL 2591 is uncertain as discussed by van 
der Tak et al. (1999) and Rygl et al. (2011). In this paper, we 
assume a distance of 1 kpc with total luminosity of 20000 Lq. 
From the radiative transfer modeling of HD0/H2^^0 with our 
toy model, we found that the azimuthal motion of the proposed 
disk-like structure is dominated by sub-Keplerian rotation which 
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is about 1/3 in strength of pure Keplerian rotation with a stellar 
mass of 16 Mq, and the radial motion is consistent with Hubble- 
law like expansion with the analogue Hubble constant of 10 km 
s"^ arcsec"^ For other distances, our results remain unchanged 
except for the relative strengths of rotation and expansion. In par- 
ticular, vro scales as 1/d and ^ scales as 1/ ^^d, so that for d = 0.5 
and 3.0 kpc, vrq is changed to 20 and 3.3 km s"^ arcsec"\ while 
^ becomes 0.5 (effective Keplerian mass = 4 Mq) and 0.2 (ef- 
fective Keplerian mass = 0.7 Mq). We suggest that the disk-like 
structure cannot undergo pure Keplerian rotation at any of these 
distances since the effective Keplerian masses are not consistent 
with the scaled total luminosities. 

6. Conclusions 

Our (sub)arcsecond resolution PdBI observations in millimeter 
continuum and in lines of HDO li,o-li,i, H2^^0 3 13-22,0 and 
SO2 12o,i2-lli,ii successfully resolve the inner thousand AU re- 
gion of the young high-mass star AFGL 2591-VLA3. We sum- 
marize our findings below. 

1. From the visibility analysis of millimeter continuum and 
molecular lines, a compact source with diameter of ~ 800 
AU at 1 kpc is found toward VLA3. Depending on the as- 
sumptions, the H2 column density, number density and mass 
are estimated to be 1 x 10^"^ cm"^, 1 x 10^ cm"^ and 0.2 Mq, 
respectively. 

2. A linear velocity gradient in the northeast-southwest direc- 
tion is found across the source, which has an offset in P.A. 
by ~ 40° to the equatorial region (P.A. 0°) defined by the 
east- west outflow. We interpret this kinematical signature as 
a combination of rotation and radial expansion in the equa- 
torial region of VLA3. 

3. By introducing a toy model with simplified parameter op- 
timization, we further constrain the rotation to be sub- 
Keplerian and the radial expansion to be Hubble-law like. 
This result is independent to the distance of the source. 

4. We speculate that sub-Keplerian rotation in the equatorial re- 
gion could be effect of magnetic field which turns the veloc- 
ity profile from pure Keplerian to sub-Keplerian. The expan- 
sion could be due to the interaction between the stellar radi- 
ation/winds and the gas in the equatorial region of VLA3. 

5. Based on the molecular line profiles, the chemical structure 
of the proposed disk appears layered, with HDO near the 
mid-plane, H2^^0 in the mid-layer and SO2 in the upper lay- 
ers, judging from their turbulent line widths. 
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